# This file estimates effects of twinning by age of the child for fathers

# install.packages("dplyr")
# install.packages("date")
# install.packages("AER")

library(dplyr)
library(date)
library(AER)

# Load data
load("mothers.rdata")

# Filter of if age of oldest is unkwown (select parents )
mothers <-
  mothers %>%
  filter(!is.na(age_oldest)) 

# split data file to each election 
# 2009 election
# Subset and create vote variable
# filter of parents with children born after the election
# loop over years 

het_results <- 
  matrix(NA, nrow = 144, ncol = 5)

for(i in 1:16){
  mothers09 <- 
    mothers %>%
    filter(! is.na(stemte_2009) & first09 > 0) %>%
    filter(age_oldest <= (14565 - (i-1)*365 + i - 1) & age_oldest > (14565 - i * 365 + i)) %>%
    mutate(vote = stemte_2009 == "Stemte",
           no_children = no_children_09)
  
  #create age dummies for mothers 
  #dichotomize twin first
  mothers09 <-
    within(mothers09, {
      age        <- ntile(date, 10)
      twin_first <- twin_first > 0
    })
  
  mothers09 <-
    mothers09 %>%
    group_by(age) %>%
    mutate(child_age = ntile(age_oldest, 2))
  
  ## Repeat for 13
  mothers13 <- 
    mothers %>%
    filter(! is.na(stemt) & first13 > 0) %>%
    filter(age_oldest <= (16028 - (i-1)*365 + i - 1) & age_oldest > (16028 - i * 365 + i)) %>%
    mutate(vote = stemt == "Stemte",
           no_children = no_children_13)
  
  mothers13 <-
    within(mothers13, {
      age        <- ntile(date, 10)
      twin_first <- twin_first > 0
    })
  
  mothers13 <-
    mothers13 %>%
    group_by(age) %>%
    mutate(child_age = ntile(age_oldest, 2))
  
  ## Repeat for 14
  mothers14 <- 
    mothers %>%
    filter(! is.na(ep_stemt) & first13 > 0) %>%
    filter(age_oldest <= (16215 -(i-1)*365 + i - 1) & age_oldest > (16215 - i * 365 + i)) %>%
    mutate(vote = ep_stemt,
           no_children = no_children_14)
  
  mothers14 <-
    within(mothers14, {
      age        <- ntile(date, 10)
      twin_first <- twin_first > 0
    })
  
  mothers14 <-
    mothers14 %>%
    group_by(age) %>%
    mutate(child_age = ntile(age_oldest, 2))
  
  # Run models 
  
  # Models with fixed effect for mothers and children's age
  mod09 <- lm(vote ~ twin_first + as.factor(age), data = mothers09)
  mod13 <- lm(vote ~ twin_first + as.factor(age), data = mothers13)
  mod14 <- lm(vote ~ twin_first + as.factor(age), data = mothers14)
  
  
  ## Instrument number of children with first child being twin
  mod09_iv  <- ivreg(vote ~ no_children + as.factor(age) | 
                       twin_first + as.factor(age), data = mothers09)
  mod13_iv  <- ivreg(vote ~ no_children + as.factor(age) | 
                       twin_first + as.factor(age), data = mothers13)
  mod14_iv  <- ivreg(vote ~ no_children + as.factor(age) | 
                       twin_first + as.factor(age), data = mothers14)
  
  ## Run first-stage models 
  mod09_fs <- lm(no_children ~ twin_first + as.factor(age), data = mothers09)
  mod13_fs <- lm(no_children ~ twin_first + as.factor(age), data = mothers13)
  mod14_fs <- lm(no_children ~ twin_first + as.factor(age), data = mothers14)
  
  # Compile results 
  
  het_results[(1 + (i-1)*9) : (i*9),] <- 
    rbind(c(summary(mod09)$coefficients[2,1:2],i,2009,1),
          c(summary(mod13)$coefficients[2,1:2],i,2013,1),
          c(summary(mod14)$coefficients[2,1:2],i,2014,1),
          c(summary(mod09_iv)$coefficients[2,1:2],i,2009,2),
          c(summary(mod13_iv)$coefficients[2,1:2],i,2013,2),
          c(summary(mod14_iv)$coefficients[2,1:2],i,2014,2),
          c(summary(mod09_fs)$coefficients[2,1:2],i,2009,3),
          c(summary(mod13_fs)$coefficients[2,1:2],i,2013,3),
          c(summary(mod14_fs)$coefficients[2,1:2],i,2014,3))
}

write.table(het_results, file = "het_results_mother.txt")
